family = "sans"),
axis.title = element_text(size = 12,
color = "black",
family = "sans"),
plot.title = element_text(face = "bold", size = 12),
strip.text = element_text(hjust = 0,
face = "bold",
size = 12,
color = "black"),
)
}
round_up <- function(x, digits = 1) {
if (isTRUE(all.equal(x, 0))) return(x)
decs <- 10^floor(log10(abs(x))) / 10^(digits - 1)
ceiling(x / decs) * decs
}
quater_rounding<-function(number){
round_up <- function(x, digits = 1) {
if (isTRUE(all.equal(x, 0))) return(x)
decs <- 10^floor(log10(abs(x))) / 10^(digits - 1)
ceiling(x / decs) * decs
}
x<-round_up(number,digits = 1)        #rounds up to 1 sigfig for order of magniturde calc and subtratiosns in if
y<-(x%>%abs()%>%log10()%>%floor())-1  #calculate order of magnatuede
digets<-number%>%                     #gets 2 and 3 diget
round_up(digits = 3)%>%
as.character()%>%
str_sub(2, 3)%>%
as.numeric()
if(digets>75){                        #if statment for quater rounding
z<-round_up(number,digits = 1)
}else if (digets>50){
z<-x-2.5*10^y
}else if (digets>25){
z<-x-5.0*10^y
}else if (digets>00){
z<-x-7.5*10^y
}else{
z<-signif(number,digits = 1)
}
z
}
half_rounding<-function(number){
round_up <- function(x, digits = 1) {
if (isTRUE(all.equal(x, 0))) return(x)
decs <- 10^floor(log10(abs(x))) / 10^(digits - 1)
ceiling(x / decs) * decs
}
x<-round_up(number,digits = 1)           #rounds up to 1 sigfig for order of magniturde calc and subtratiosns in if
y<-(x%>%abs()%>%log10()%>%floor())-1     #calculate order of magnatuede
digets<-number%>%                        #gets 2 and 3 diget
round_up(digits = 3)%>%
as.character()%>%
str_sub(2, 3)%>%
as.numeric()
if(digets>50){                          #if statment for half rounding
z<-round_up(number,digits = 1)
}else if (digets>00){
z<-x-5.0*10^y
}else{
z<-signif(number,digits = 1)
}
z
}
graph_rounding<-function(number){
round_up <- function(x, digits = 1) {
if (isTRUE(all.equal(x, 0))) return(x)
decs <- 10^floor(log10(abs(x))) / 10^(digits - 1)
ceiling(x / decs) * decs
}
x<-round_up(number,digits = 1)               #rounds up to 1 sigfig for order of magniturde calc and subtratiosns in if
y<-(x%>%abs()%>%log10()%>%floor())-1         #calculate order of magnatuede
digets<-number%>%                            #exstracts 2 and 3 diget
round_up(digits = 3)%>%
as.character()%>%
str_sub(2, 3)%>%
as.numeric() #
diget_1<-number%>%                           #exstracts 1st diget
round_up(digits = 3)%>%
as.character()%>%
str_sub(1,1)%>%
as.numeric()
if(digets>50 | diget_1>=4){                 #if statment for rounding to nice graphs
z<-round_up(number,digits = 1)
}else if (digets>00){
z<-x-5.0*10^y
}else{
z<-signif(number,digits = 1)
}
z
}
limmet_maker<-function(a=0,b=p1){
graph_rounding<-function(number){
round_up <- function(x, digits = 1) {
if (isTRUE(all.equal(x, 0))) return(x)
decs <- 10^floor(log10(abs(x))) / 10^(digits - 1)
ceiling(x / decs) * decs
}
x<-round_up(number,digits = 1)               #rounds up to 1 sigfig for order of magniturde calc and subtratiosns in if
y<-(x%>%abs()%>%log10()%>%floor())-1         #calculate order of magnatuede
digets<-number%>%                            #exstracts 2 and 3 diget
round_up(digits = 3)%>%
as.character()%>%
str_sub(2, 3)%>%
as.numeric() #
diget_1<-number%>%                           #exstracts 1st diget
round_up(digits = 3)%>%
as.character()%>%
str_sub(1,1)%>%
as.numeric()
if(digets>50 | diget_1>=4){                 #if statment for rounding to nice graphs
z<-round_up(number,digits = 1)
}else if (digets>00){
z<-x-5.0*10^y
}else{
z<-signif(number,digits = 1)
}
z
}
c(a,                                      #sets lower limmet
layer_data(b)%>%                        #retreaves data from plot
filter(!is.na(x) & !is.na(y))%>%      #filters to ensure no duplicate values are present
summarise(max(x,y,na.rm = TRUE))%>%   #gets max value in both x and y aes
pull()%>%                             #pulls value out of tibble
graph_rounding())                     #rounds to the apropreat value for nice looking graphs.
}
#slope metric construction
#see https://www.tidymodels.org/learn/develop/metrics/ to understand how to construct your own metric.
#the below code produces the metric slope for metric_set
#dependancy rlang and yardstick, dplyr
slope_impl<-function(truth,estimate,case_weights = NULL){
lm(truth ~ estimate)%>%
tidy()%>%
filter(term!="(Intercept)")%>%
pull(estimate)
}
slope_vec <- function(truth, estimate, na_rm = TRUE, case_weights = NULL, ...) {
check_numeric_metric(truth, estimate, case_weights)
if (na_rm) {
result <- yardstick_remove_missing(truth, estimate, case_weights)
truth <- result$truth
estimate <- result$estimate
case_weights <- result$case_weights
} else if (yardstick_any_missing(truth, estimate, case_weights)) {
return(NA_real_)
}
slope_impl(truth, estimate, case_weights = case_weights)
}
slope <- function(data, ...) {
UseMethod("slope")
}
slope <- new_numeric_metric(slope, direction = "minimize")
slope.data.frame <- function(data, truth, estimate, na_rm = TRUE, case_weights = NULL, ...) {
numeric_metric_summarizer(
name = "slope",
fn = slope_vec,
data = data,
truth = !!enquo(truth),
estimate = !!enquo(estimate),
na_rm = na_rm,
case_weights = !!enquo(case_weights)
)
}
#intercept metric construction
#see https://www.tidymodels.org/learn/develop/metrics/ to understand how to construct your own metric.
#the below code produces the metric intercept for metric_set
#dependancy rlang and yardstick, dplyr
intercept_impl<-function(truth,estimate,case_weights = NULL){
lm(truth ~ estimate)%>%
tidy()%>%
filter(term=="(Intercept)")%>%
pull(estimate)
}
intercept_vec <- function(truth, estimate, na_rm = TRUE, case_weights = NULL, ...) {
check_numeric_metric(truth, estimate, case_weights)
if (na_rm) {
result <- yardstick_remove_missing(truth, estimate, case_weights)
truth <- result$truth
estimate <- result$estimate
case_weights <- result$case_weights
} else if (yardstick_any_missing(truth, estimate, case_weights)) {
return(NA_real_)
}
intercept_impl(truth, estimate, case_weights = case_weights)
}
intercept <- function(data, ...) {
UseMethod("intercept")
}
intercept <- new_numeric_metric(intercept, direction = "minimize")
intercept.data.frame <- function(data, truth, estimate, na_rm = TRUE, case_weights = NULL, ...) {
numeric_metric_summarizer(
name = "intercept",
fn = intercept_vec,
data = data,
truth = !!enquo(truth),
estimate = !!enquo(estimate),
na_rm = na_rm,
case_weights = !!enquo(case_weights)
)
}
setwd("C:/Users/s299038/OneDrive - Cranfield University/1. Rapid Measurement Tools - Niall/Chapter 3 - XRF Analysis Accuracy on Slag/Data Analysis/combining for analysis")
pivoted_data<-read_csv("filterd_data.csv")
setwd("C:/Users/s299038/OneDrive - Cranfield University/1. Rapid Measurement Tools - Niall/Chapter 3 - XRF Analysis Accuracy on Slag/Data Analysis/Fig_J wet sample vs dry sample")
setwd("C:/Users/s299038/OneDrive - Cranfield University/1. Rapid Measurement Tools - Niall/Chapter 3 - XRF Analysis Accuracy on Slag/Data Analysis/Apendix")
element_orders<-read_csv("instrument_elements.csv")
print(element_orders,n=26)
elements_conc<-element_orders$Conc_Order_reverse%>%na.omit()
elements_atomic<-element_orders$Atomic_Order%>%na.omit()
elements_alphabetical<-element_orders$Alphabetic_Order%>%na.omit()
setwd("C:/Users/s299038/OneDrive - Cranfield University/1. Rapid Measurement Tools - Niall/Chapter 3 - XRF Analysis Accuracy on Slag/Data Analysis/Fig_J wet sample vs dry sample")
dried_data<-pivoted_data%>%
filter(Prosessing_Step=="dried")%>%
select(ID,element,mean,SD,RSD)%>%
rename(mean_dried=mean,SD_dried=SD,RSD_dried=RSD)
wet_data<-pivoted_data%>%
filter(Prosessing_Step=="sived")%>%
select(!Prosessing_Step)%>%
rename(mean_wet=mean,SD_wet=SD,RSD_wet=RSD)
grind_data<-pivoted_data%>%
filter(Prosessing_Step=="ground")%>%
select(!Prosessing_Step)%>%
rename(mean_grind=mean,SD_grind=SD,RSD_grind=RSD)
ign_data<-pivoted_data%>%
filter(Prosessing_Step=="ignited")%>%
select(!Prosessing_Step)%>%
rename(mean_ign=mean,SD_ign=SD,RSD_ign=RSD)
joined_data<-
full_join(wet_data,dried_data,by=c("ID","element"))%>%
full_join(grind_data)%>%
full_join(ign_data)%>%
mutate(element = fct_relevel(element,elements_conc),
Depth = as.factor(Depth),
matrix = as.factor(matrix))
#obs_pred scales
element<-levels(joined_data$element)
color_vector <- brewer.pal(length(element), "Paired")
for (i in 1:length(element)) {
e<-element[i]
f<-color_vector[i]
p1<-joined_data%>%
filter(element==e)%>%
ggplot(aes(x=mean_dried,y=mean_wet))
p1<-p1+ #this break in the ggplot is needed so limmet_maker can acsess the plot data in p1
geom_smooth(method=lm,size=1,color = "blue",alpha = 1,se = FALSE)+
geom_errorbar(aes(xmin = mean_dried - SD_dried, xmax = mean_dried + SD_dried),
width = 0.1,
color = "black",
alpha = .5)+
geom_errorbar(aes(ymin = mean_wet - SD_wet, ymax = mean_wet + SD_wet,),
width = 0.1,
color = "black",
alpha = .5)+
geom_point(size=1,aes(color = Depth,shape=matrix))+
geom_abline(aes(intercept=0, slope=1),linetype="dashed")+
coord_cartesian(xlim = limmet_maker(), ylim=limmet_maker())+
scale_y_continuous_excel(name =NULL)+
scale_x_continuous_excel(name =NULL)+
labs(title = e)+
my_theme()+
theme(
axis.text.x = element_text(angle=25,hjust = 1,)
)
# Add the dataframe to the list
if(i==1){
plots<-list()
}
plots[[e]] <- p1
if(i==1){
patchwork_plot<-p1
}else if(i==length(element)){
patchwork_plot<-patchwork_plot+p1
rm(p1)
rm(e)
rm(element)
rm(i)
rm(color_vector)
rm(f)
}else{
patchwork_plot<-patchwork_plot+p1
}
}
patchwork_plot+plot_layout(ncol = 3)
ggsave("Fig_J_moisture_faceted_errorbars.png",width = 160,height = (160/3)*4, units = "mm")#A4 margens square use hight to calc for grided plots
ggsave("Fig_J_moisture_faceted_errorbars.svg",width = 160,height = (160/3)*4, units = "mm")#A4 margens square use hight to calc for grided plots
#obs_pred scales
element<-levels(joined_data$element)
color_vector <- brewer.pal(length(element), "Paired")
for (i in 1:length(element)) {
e<-element[i]
f<-color_vector[i]
p1<-joined_data%>%
filter(element==e)%>%
ggplot(aes(x=mean_ign,y=mean_grind))
p1<-p1+ #this break in the ggplot is needed so limmet_maker can acsess the plot data in p1
geom_smooth(method=lm,size=1,color = "black",alpha = 1,se = FALSE)+
geom_errorbar(aes(xmin = mean_ign - SD_ign, xmax = mean_ign + SD_ign),
width = 0.1,
color = "black",
alpha = .5)+
geom_errorbar(aes(ymin = mean_grind - SD_grind, ymax = mean_grind + SD_grind,),
width = 0.1,
color = "black",
alpha = .5)+
geom_point(size=1,aes(color = Depth,shape=matrix))+
geom_abline(aes(intercept=0, slope=1),linetype="dashed")+
coord_cartesian(xlim = limmet_maker(), ylim=limmet_maker())+
scale_y_continuous_excel(name =NULL)+
scale_x_continuous(name =NULL)+
labs(title = e)+
my_theme()+
theme(
axis.text.x = element_text(angle=25,hjust = 1,)
)
# Add the dataframe to the list
if(i==1){
plots<-list()
}
plots[[e]] <- p1
if(i==1){
patchwork_plot<-p1
}else if(i==length(element)){
patchwork_plot<-patchwork_plot+p1
rm(p1)
rm(e)
rm(element)
rm(i)
rm(color_vector)
rm(f)
}else{
patchwork_plot<-patchwork_plot+p1
}
}
plots[[4]]
patchwork_plot+plot_layout(ncol = 3)
ggsave("Fig_J_organic_faceted_errorbars.png",
device = png,
width = 160,
height = (160/3)*4,
units = "mm")#A4 margens square use hight to calc for grided plots
ggsave("Fig_J_organic_faceted_errorbars.svg",
device = svg,
width = 160,
height = (160/3)*4,
units = "mm")#A4 margens square use hight to calc for grided plots
#obs_pred scales
element<-levels(joined_data$element)
color_vector <- brewer.pal(length(element), "Paired")
for (i in 1:length(element)) {
e<-element[i]
f<-color_vector[i]
p1<-joined_data%>%
filter(element==e)%>%
ggplot(aes(x=mean_ign,y=mean_grind))
p1<-p1+ #this break in the ggplot is needed so limmet_maker can acsess the plot data in p1
geom_smooth(method=lm,size=1,color = "black",alpha = 1,se = FALSE)+
geom_errorbar(aes(xmin = mean_ign - SD_ign, xmax = mean_ign + SD_ign),
width = 0.1,
color = "black",
alpha = .5)+
geom_errorbar(aes(ymin = mean_grind - SD_grind, ymax = mean_grind + SD_grind,),
width = 0.1,
color = "black",
alpha = .5)+
geom_point(size=1,aes(color = Depth,shape=matrix))+
geom_abline(aes(intercept=0, slope=1),linetype="dashed")+
coord_cartesian(xlim = limmet_maker(), ylim=limmet_maker())+
scale_y_continuous_excel(name =NULL)+
scale_x_continuous_excel(name =NULL)+
labs(title = e)+
my_theme()+
theme(
axis.text.x = element_text(angle=25,hjust = 1,)
)
# Add the dataframe to the list
if(i==1){
plots<-list()
}
plots[[e]] <- p1
if(i==1){
patchwork_plot<-p1
}else if(i==length(element)){
patchwork_plot<-patchwork_plot+p1
rm(p1)
rm(e)
rm(element)
rm(i)
rm(color_vector)
rm(f)
}else{
patchwork_plot<-patchwork_plot+p1
}
}
plots[[4]]
patchwork_plot+plot_layout(ncol = 3)
ggsave("Fig_J_organic_faceted_errorbars.png",
device = png,
width = 160,
height = (160/3)*4,
units = "mm")#A4 margens square use hight to calc for grided plots
ggsave("Fig_J_organic_faceted_errorbars.svg",
device = svg,
width = 160,
height = (160/3)*4,
units = "mm")#A4 margens square use hight to calc for grided plots
sumerising<-function(data){
#set number of levels needed for loop
element_nleval<-nlevels(as.factor(data$element))
element_levals<-levels(as.factor(data$element))
for(l in 1:element_nleval){
#pick the eliment of this loop
eliment<-element_levals[l]
#filter data
filterd_data<-data%>%filter(element==eliment)
#creat linnear model for slope and intersept
linnear_model<-lm(mean_wet ~ mean_dried,data = filterd_data)
#exstract slop/intersept and wrangel in to tibble
slope_intersept<-bind_rows(linnear_model[[1]])%>%
rename(slope=mean_dried,intercept=`(Intercept)`)
#exstract confidence intervals and wrangel in to tibbel
c_1<-confint(linnear_model, level = 0.95)%>%
as_tibble()%>%
slice(1)%>%
mutate(across(everything(), ~ round(.,digits = 0)))%>%
summarise("Intersept 95% CI"=str_c("(",`2.5 %`,", ",`97.5 %`,")"))
c_2<-confint(linnear_model, level = 0.95)%>%
as_tibble()%>%
slice(2)%>%
mutate(across(everything(), ~ round(.,digits = 3)))%>%
summarise("Slope 95% CI"=str_c("(",`2.5 %`,", ",`97.5 %`,")"))
CI<-cbind(c_1,c_2)
#calculate n
N<-filterd_data%>%summarise(n=n())
#build metric set
metrics<-metric_set(rpd,rsq,rmse,rpiq,ccc,mae,rsq_trad)
#calculat meterics and organise and add slope/intersept and add names and move to front
element_stats<-filterd_data%>%
metrics(truth = mean_dried,estimate = mean_wet)%>%
pivot_wider(id_cols = .estimator,names_from = .metric,values_from = .estimate)%>%
mutate(N,
slope_intersept,
CI,
element=eliment)%>%
select(!.estimator)%>%
relocate(element)
#if function to combine data in to 1 df
if(c(l==1)){
sum_stats<-bind_rows(element_stats)
}else{
sum_stats<-bind_rows(sum_stats,element_stats)
}
#close loop
}
#print data
sum_stats
}
wet_vs_dry_tabel<-sumerising(joined_data)
View(wet_vs_dry_tabel)
wet_vs_dry_tabel%>%write_csv("Fig_J_moisture_suplimentry_tabel.csv")
sumerising<-function(data){
#set number of levels needed for loop
element_nleval<-nlevels(as.factor(data$element))
element_levals<-levels(as.factor(data$element))
for(l in 1:element_nleval){
#pick the eliment of this loop
eliment<-element_levals[l]
#filter data
filterd_data<-data%>%filter(element==eliment)
#creat linnear model for slope and intersept
linnear_model<-lm(mean_grind ~ mean_ign,data = filterd_data)
#exstract slop/intersept and wrangel in to tibble
slope_intersept<-bind_rows(linnear_model[[1]])%>%
rename(slope=mean_ign,intercept=`(Intercept)`)
#exstract confidence intervals and wrangel in to tibbel
c_1<-confint(linnear_model, level = 0.95)%>%
as_tibble()%>%
slice(1)%>%
mutate(across(everything(), ~ round(.,digits = 0)))%>%
summarise("Intersept 95% CI"=str_c("(",`2.5 %`,", ",`97.5 %`,")"))
c_2<-confint(linnear_model, level = 0.95)%>%
as_tibble()%>%
slice(2)%>%
mutate(across(everything(), ~ round(.,digits = 3)))%>%
summarise("Slope 95% CI"=str_c("(",`2.5 %`,", ",`97.5 %`,")"))
CI<-cbind(c_1,c_2)
#calculate n
N<-filterd_data%>%summarise(n=n())
#build metric set
metrics<-metric_set(rpd,rsq,rmse,rpiq,ccc,mae,rsq_trad)
#calculat meterics and organise and add slope/intersept and add names and move to front
element_stats<-filterd_data%>%
metrics(truth = mean_ign,estimate = mean_grind)%>%
pivot_wider(id_cols = .estimator,names_from = .metric,values_from = .estimate)%>%
mutate(N,
slope_intersept,
CI,
element=eliment)%>%
select(!.estimator)%>%
relocate(element)
#if function to combine data in to 1 df
if(c(l==1)){
sum_stats<-bind_rows(element_stats)
}else{
sum_stats<-bind_rows(sum_stats,element_stats)
}
#close loop
}
#print data
sum_stats
}
ground_vs_ignited_tabel<-sumerising(joined_data)
View(ground_vs_ignited_tabel)
ground_vs_ignited_tabel%>%write_csv("Fig_J_organic_suplimentry_tabel.csv")
